The existence of relativistic stars in f(R) gravity 



o 
o 

(N 



Amol Upadhye 1 and Wayne Hu 1,2 

1 Kavli Institute for Cosmological Physics, Enrico Fermi Institute, University of Chicago, Chicago, IL 60637 
^Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637 

(Dated: May 25, 2009) 

We refute recent claims in the literature that stars with relativistically deep potentials cannot 
exist in f(R) gravity. Numerical examples of stable stars, including relativistic (Gh'h/r* ~ 0.1), 
constant density stars, are studied. As a star is made larger, non-linear "chameleon" effects screen 
much of the star's mass, stabilizing gravity at the stellar center. Furthermore, we show that the 
onset of this chameleon screening is unrelated to strong gravity. At large central pressures P > p/3, 
f(R) gravity, like general relativity, does have a maximum gravitational potential, but at a slightly 
smaller value: GM*/r*j max = 0.345 < 4/9 for constant density and one choice of parameters. This 
difference is associated with negative central curvature R under general relativity not being accessed 
in the f(R) model, but does not apply to any known astrophysical object. 
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I. INTRODUCTION 

The discovery that the expansion of the universe is ac- 
celerating spurred a search for theoretical models which 
could explain this phenomenon. The simplest explana- 
tion, the cosmological constant, requires extreme fine 
tuning in order to explain its smallncss as well as its 
closeness to today's matter density. This motivates the 
search for alternative explanations for the cosmic accel- 
eration. These alternatives fall into two broad classes. 
In the first class, a new field, known as a "dark energy" , 
comes to dominate the universe at recent times, prevent- 
ing the Hubble parameter from falling as rapidly as it 
would in a matter-dominated universe. The second class 
alters gravity itself, with modifications on large scales 
causing the universe to deviate from its expected decel- 
eration today. 

Modified gravity explanations are highly constrained 
by our knowledge of gravity on small scales. These hur- 
dles include equivalence principle tests [l|, 0] , solar sys- 
tem measurements 0] , and the stability of gravitationally 
bound systems such as stars [3, We consider f(R) 
gravity, a theory in which the Ricci scalar R is replaced 
by some function f(R) in the action for gravity [a0,[l]. 
While it has been shown that some models are consis- 
tent with solar system measurements @ and the sta- 
bility of non-relativistic stars ||, Kobayashi and Maeda 
(hereafter KM) have recently claimed that relativis- 
tic stars are unstable in related f(R) models These 
arguments point to the existence of a curvature singular- 
ity in cosmologically viable theories; as the scalar field 
<j> = df/dR — ► 1, R — ► oo. They claim that relativistic 
stars, with GM+fr* ~ 0.1, push <j) right into the curva- 
ture singularity, meaning that stars which we know to 
exist could not in f(R) gravity. Other works, while not 
directly disproving this claim, have argued that the sin- 
gularity may be avoided by choosing a different equation 
of state [12|, or a model in which a divergence in the 
scalar field potential ensures that <f> avoids the singular- 
ity (HQ. 

Here we show, through numerical computation as well 



as analytical argument that highly relativistic stars do in- 
deed exist in f(R) gravity. Existence does not hinge on 
a specific equation of state or choice of f(R) but rather 
the non-linearity of the field equations. The onset of 
non-linearity causes the field to stop changing with the 
potential via the so-called chameleon effect. Thereafter 
deviations in <f> from its background value are determined 
only by a small portion of the stellar mass and the cur- 
vature singularity is never reached in a static star. Non- 
linearity in the equations of motion make the numeri- 
cal solutions difficult to attain, which has obscured these 
points in the literature. Nonetheless, we have numerically 
confirmed the existence of ultra-relativistic stars with po- 
tentials GM+/r+ up to 0.345 and central pressures much 
greater than their energy densities. 

Furthermore, we show that the onset of non-linear 
chameleon effects has nothing to do with strong grav- 
ity. They will generically arise when the gravitational 
potential GM+/r+ exceeds the field distance between the 
background value and the curvature singularity, which 
depends on the f(R) function itself. In KM [1(|, this 
distance was taken to be of order 0.1. In fact this dis- 
tance must be < 10~ 6 — 10~ 5 to remain compatible with 
local tests of gravity due to the finite extent of our galaxy 
1- 

The paper is organized as follows. After introducing 
f(R) theory and its application to stars in Sec. [Til we 
present our numerical solutions in Sec. IIII1 including rel- 
ativistic as well as non-relativistic stars. Sec. lIVI cmplovs 
analytic arguments in the linear and non-linear regimes 
to elucidate how chameleon screening allows the field to 
avoid the curvature singularity. We conclude in Sec. |Vj 



II. FORMALISM 

We briefly review the equations governing f(R) theory 
in Sec. Ill Al apply them to stellar, static, spherically sym- 
metric cases in Sec. lIIBl and specialize to the Starobinsky 
f(R) model [n| in Sec.HTCl 
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A. f(R) theory 



the equation of hydrostatic equilibrium, 



Replacing the Ricci scalar R in the Einstein-Hilbert ac- 
tion defining general relativity by a function f(R) results 
in the action 



S 



(W) 



(1) 



Vf67rG 

and the modified Einstein equation 

4>R„u - + g^Ucj) - -fg^ = 8nGT^ v . (2) 

The quantity = df jdR behaves as a scalar field coupled 
to matter and the metric, as can be seen by taking the 
trace of Eq. (0), 



^T + i[2/W))-W)] 



8ttG, 
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d<t> 



eft 



(3) 



where R is now an implicit function of </>. Evidently the 
scalar is a chameleon field [HI, ; its self interaction 

V(4>) and its coupling to matter give the field a constant 
value in a medium of constant T = = —p + 3P which 
also determines its mass. We will see that non-linear ef- 
fects associated with changes in T between two differ- 
ent media are crucial for understanding f(R) solutions 
in stars. 



B. Spherical stars 

Since we intend to study stars, we assume a spherically 
symmetric metric, 



ds A 



-N{r)dt 2 



dr 2 
W. 



+ r 2 (d0 2 + sin 2 9dtp 2 ). (4) 



With this metric, the field equation ^ for static solu- 
tions becomes 



r 2JV 25 



(5) 

The system is completed by the (tt) and (rr) components 
of the modified Einstein equations, 
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-torGp-l<l>R(<l>) + lf(R(<j>)), (6) 



-1 + -B 
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2 - + ^)<t>'B 
r 2N 



MGP-^R{4>) + \f{R{4>)), (7) 



P> = - m ip + P), 



(8) 



and equation of state p = p(P) for the matter. We fol- 
low KM and assume a constant density, p(P) = p c , 
and central pressure P c , but this may easily be gener- 
alized. Here, and throughout the paper, primes denote 
derivatives with respect to r. 

Boundary conditions for <//, P, N, and B can be spec- 
ified at the center. Continuity of the gradient of 4> at 
the center of the star requires <fi'(0). P{0) is set to a 
specified central pressure P c . In order to facilitate com- 
parison with KM, we take N(0) = B(0) = 1, amounting 
to a rescaling of the time coordinate. 

The remaining boundary condition for the field is more 
complicated. We again follow KM and take the exte- 
rior of the star to be empty save for the <p field and 
any effective cosmological constant that its value implies. 
In general relativity, the exterior metric would be the 
Schwarzschild-de Sitter spacetime, which has a horizon 
J"n = J"b ra \/3/A where N(r^) = and B(rs) = 0. 
In f(R) gravity, N and B do not necessarily vanish 
at the same position. However, the curvature invariant 
Rap~tsR a ^ S diverges at if 7^ vq . In order to ensure 
that a solution of the equations of motion yields a well- 
behaved cosmology at the horizon, our final boundary 
condition must be tn = ^b- The field equation (J5]) and 
linearity in N and B around tn indicates that this bound- 
ary condition is equivalent to taking B'<fi' — V^| r - >?- N = 0. 

For a star with specified p c and P c , we find <j>(r), N(r), 
B{r), and P{r) using a shooting method. We begin by 
guessing a central value for the field, (f>(0) = <p c . With 
the boundary and normalization conditions above, we 
can solve the above system of equations, thereby find- 
ing B'4>' — Vif, at the horizon r^. Suppose we have two 
guesses, <p c - and <f> c ,+, for which B'<j>' — V t< j, < and > 0, 
respectively. Since the system of equations ([5][5]) provides 
a continuous mapping from <p c to B' <p' — at the hori- 
zon, there must be a value <f> c between C ._ and <p c .+ f° r 
which the boundary condition at the horizon is satisfied. 
We find an improved estimate of the true solution by 
iteratively bisecting this interval until the dimcnsionlcss 
boundary condition, b((j) c ) = r^ifi — V t< f,/B'), evaluated 
when TV drops to 10~ 15 , is within 1CP 5 of zero. Note that 
for stars that are much smaller than the horizon, <fi ss (jj^s 
and nearly constant, so <p' and both vanish at tn- 

Even in cases where it is numerically difficult to find 
the solution to the required precision, a solution should 
still exist if both under and overshoot solutions also ex- 
ist. This is simply the overshoot-undershoot argument of 
[lH , [l9j . applied to spherical stars in an f(R) analogue 
of a Schwarzschild-de Sitter background, rather than to 
spherical bubbles in a de Sitter background. Similarly 
to those references, our field equation ([5]) looks like the 
equation of motion of a particle at "position" 0, as a 
function of "time" r, in a potential U(<j>) = —V(<p), with 
a time-dependent friction term ~ 2/r, driven by a time- 
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FIG. 1: Inverted potential U(</>) = —V(<f>), up to an arbitrary 
additive constant, for n = 1 and xi = 3.6. 



dependent force T — 8irGT/3. For stars much smaller 
than the horizon, the vanishing of dV/d<p at the hori- 
zon corresponds to the field stopping at the maximum of 
U{4>), corresponding to the false vacuum in [l9| . 
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FIG. 2: Properties of the potential U((fi) as a function of 
x\ = Rds/Ro for n — 1. Increasing xi — that is, decreasing Rq 
at fixed Rds = 4A b s — pushes the maximum of the potential 
closer to the curvature singularity. Note that the effective 
mass mas vanishes as xi — > v^3- 



2nA/(l - cf>), so 



C. Starobinsky /(i?) model 

The above discussion applies to all f(R) models. In or- 
der to proceed, we must specify a model. To test the KM 
claim, we follow them in choosing the model of Starobin- 

sky pa, 



f(R) = R + XRq 



1 R2 
n 



(9) 



where n, A, and i?o are free parameters, and the field 



1 - 2n\ 
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by 



The potential U(4>), defined up to an additive constant 



dU 



1 



-R 



2XR Q 



1 
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l)R*/Rl 



(1 + RVR 2 ) 



,(H) 



is shown in Fig. Q] Note that <f> — 1 corresponds to a 
curvature singularity, R = oo. This is a common feature 
of /(i?) where the modification to the Einstein-Hilbcrt 
action vanishes at high curvature, including Starobinsky 
models with n > and the broken power law models of 
Hu and Sawicki For </> near the singularity, R 2> Ro, 
so U'((f>) w -R((j))/B. Eq. (HHD implies (R/R ) 2n+1 « 



-3^0 



2nA 



The effective chameleon mass 
singularity is 



(12) 



-E/"(0) near the 



d(j)/dR 
Ro 



R 



Rq 



6n(2n+ 1)A ^1- 



277.A 



6n(2n + 1)A 

2n + l 



i?o 



2n+2 



(13) 



As the singularity is approached, the effective mass and 
the slope of the potential diverge. 

Also, U has a maximum at some field value 0ds- Far 
from a star, <f> will approach c/>ds, corresponding to a 
background, de Sitter-like spacctime. In order to make 
the de Sitter background universe resemble ours, we re- 
quire that Rds = R{<j>ds) = 4A Q bs, where A D b s is the ob- 
served value of the cosmological constant. With this con- 
straint, the two constants A and i?o, as well as the posi- 
tion of the maximum c/>ds , are specified by the parameter 
x\ = Rds/R-o, as shown in Fig. O The figure also shows 
the effective chameleon mass at 4>ds, "^ds = — ^"(<Ms)- 

Note that KM chose a value of x\ = 0(1) for their 
tests and hence a value of 1 — </><js ~ 0.1. We shall see 
that that choice is responsible for the appearance of non- 
linear effects only in relativistic stars, but is not viable 
due to solar system tests of gravity 0]. 
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FIG. 3: Chameleon field <j> (left) and its derivative (j)' (right) in a non-relativistic star with P c /p c — 0.01 and GM*/r+ = 0.0139, 
as well as a relativistic star with P c /p c = 0.1 and GM 4 /r t = 0.131. We have assumed n = 1, x\ = 3.6, and p c = 100pA- 
The solid vertical line denotes the position of the horizon, and the dotted lines are the stellar radii, with the larger radius 
corresponding to the relativistic star. In the plot of <f>' , on the right, thick lines represent positive 4>' an d thin lines represent 
negative <j)' . H/y = a/A/3 ~ 1/Vn is the horizon scale in the analogous general relativistic spacetime. 



III. NUMERICAL SOLUTIONS 



center of the star, we have 



We discuss numerical considerations in Sec. IIII Al that 
lead to the choice of stellar parameters for which we give 
numerical solutions in Sec. IIII Bl These examples are cho- 
sen to have gravitational potentials that are comparable 
to or exceed astrophysical stars albeit at a much lower 
density. These considerations and those of the next sec- 
tion suggest that stars with realistic densities should also 
exist, even though these cases are numerically intractable 
with our techniques. 



A. Numerical considerations 

A chameleon field can change on distance scales of or- 
der its Compton wavelength m^ 1 (sec Eq. IT3"|). Because 
Yukawa-like error modes grow as e m, *' r /r, numerical so- 
lutions become rapidly intractable as this scale becomes 
much smaller than the computational domain. 

In order to study the chameleon field numerically in a 
star of radius r+ in a brute force implementation, we must 
use more than m0r max integration steps. For our calcula- 
tion r max = r]y, the horizon, but in general r max 3> r+ in 
order to match an exterior solution. For small stars, the 
chameleon field is in the linear regime, where 4> is only 
slightly perturbed from </>ds, and m^(^) ~ mas [20L l2lf. 
Fig. [2] shows that m^s lies within a few orders of magni- 
tude of the horizon scale H\ = \/ A/3 for a large range 
of x\. Because Hpj"* "C 1, the chameleon equations are 
numerically tractable for small stars. 

On the other hand, for large stars <f> will be in the non- 
linear regime, R ^> Ro. Assuming R « 8irG(p c — 3P C ), 
for which <p c minimizes the effective potential V s at the 
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/2 f SnGpt 
Ro 

D \ 1/2 /o \ n+l/2 

Pc) V Ro 



1/2 



(14) 



where the non-relativistic approximations P c <C p c and 
r* 2 = 12P c /(87rGp 2 ) have been used. For x\ of order 
unity, Ro ~ 8nGp\. Thus, if we want m^r* < 1000, 
then we must have p c < 100/9A for n = 1, and even lower 
p c for larger n. 

On top of these issues, the shooting technique exacer- 
bates the difficulty in finding solutions that satisfy the ex- 
terior boundary condition. Again because of the Yukawa- 
like nature of the solutions, the central field value </> c 
must be exponentially tuned to give the correct bound- 
ary value (see Sec. HVBI for a more extended discussion). 
For example, consider the star with n = 1, x\ = 3.6, 
p c = 100pA, and P c /p c = 0.1, with the boundary condi- 
tion b{(j) c ) = r^{4>' — Vfi/B') = at the horizon. We find 
numerically that, if we want |6(</> c )| < 10 -5 at the hori- 
zon, then <p c must be tuned to within 1.5 x 10~ 37 of its 
correct value. Requiring |6((/> c )| < 10 -10 means tuning 
4> c to within 1.5 x 10 -42 of its correct value. To avoid 
these issues, relaxation methods can be applied instead 
@ but we choose a shooting method to test the KM claim 
directly. 

KM note that denser stars are more numerically dif- 
ficult, and argues that only the gravitational potential 
is relevant as a measure of the star's size; an instability 
in a large, low-density star should persist in a smaller, 
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denser star of the same gravitational potential. KM then 
go on to choose p c = 2 x 10 6 pA, f ar from the density of 
a typical star. However, even this density is too large. 
By Eq. p4[) . integration of the equations of motion for 
a star with this density would require billions of inte- 
gration steps for a brute force approach, and shooting 
compounds this problem by requiring a large number of 
these solutions to iterate to the proper boundary condi- 
tions. 

This numerical difficulty appears as the field fluctua- 
tions transition to the non-linear regime where the Comp- 
ton wavelength shrinks substantially in the stellar inte- 
rior. If we choose a large, constant density and gradually 
increase the "size" GM+/r+ of the star, then the integra- 
tion of the equations of motion will rapidly become more 
difficult as the non-linear regime is approached. In the 
case of stars needing billions of integration steps, trunca- 
tion errors may make integration impossible. Nonetheless 
this numerical difficulty does not imply that solutions do 
not exist. 



B. Example Solutions 

Contrary to the claim of KM [13], we find that 
stellar solutions exist at large gravitational potentials, 
GM+jr* > 0.1. As discussed above, we set n — 1 and 
p c = 100/3A in order to keep the problem tractable. We 
also set x\ = 3.6, unless otherwise specified, for ease of 
comparison with KM. In each case we have integrated 
the equations of motion (JSJH]) directly, using a Rungc- 
Kutta-Fehlberg mixed 4th/5th order algorithm with vari- 
able step sizes. We also emp loy the arbitrary precision 
arithmetic package CLN [22|], and all our computations 
use at least 50 decimal places. Fig. [3] shows cf>(r) and 
4>'(r) for a non-relativistic star with P c = 10 _2 p c and 
GM 4 /r t = 0.0139, as well as a relativistic star with 
P c = 10~Vc and GM*/r* = 0.131. 

For comparison purposes, we have also attempted to 
extend our computations to higher a central density 
p c = 2 x 10 6 pa, as in KM. For unsaturated stars with 
pressures of 10 _4 p c and 0.05p c , we find solutions that 
agree with those in KM. Stars in the non-linear regime 
are numerically intractable, as expected. For P c = 0.1/? c , 
numerical instabilities prevent us from following 4>{r) be- 
yond r ~ 10 -5 7v Our <fi'(r) oscillates with an increasing 
amplitude about some central value, and eventually flies 
off to large positive or negative values, leading to an un- 
dershoot or an overshoot, respectively 

Our lower central density choice makes numerical so- 
lutions tractable for the full range of central pressures. 
It is evident from Fig. [3] (left) that increasing P c from 
10~ 2 p c to lO^Vc (with corresponding increases in the 
gravitational potential) causes a greater perturbation in 
the field <f> c from its background value toward the curva- 
ture singularity. 

The mechanical analogy between (f> and a particle in 
the potential U((f>) is helpful in understanding the behav- 



ior of the field and why it does not hit the curvature sin- 
gularity as P c is raised further. From the shape of the po- 
tential shown in Fig.[H we see that larger (j) c corresponds 
to a larger magnitude of the slope U $ of the potential, 
and Eq. (fT"2| implies that this slope diverges as <f) c — ► 1. 
On the other hand, a larger pressure corresponds to a 
smaller value of the "force" term = 8irG(p — 3P)/3 
in the field equation of motion T pushes the field to 
lower values, allowing it to roll towards the peak of 
the potential. 

Since increasing P c causes the slope \U^\ to rise and 
the force \J-\ to fall, there must be some central pressure 
beyond which the slope of the potential overwhelms the 
force, and the field cannot decrease near the center. From 
Eq. ([5]) we see that the threshold value (j>t at which the 
force at r = precisely cancels the slope of the potential 
is the minimum of the effective potential V c s, 







dVcs 



dcj> 

t « 1 - 2n\ 



dU 

~d4> 



8ttG 



(p c - 3P C ) 



Po 



8ttG( P c - 3P C ) 



(15) 



One argument against the existence of relativistic stars 
is that it is possible to increase P c until <j) c > 4>t- In 
such a star, 4> will increase with r near the stellar cen- 
ter, approaching the singularity 0=1. However, the 
curve corresponding to P c / p c = 0.1 in Fig. [3] (right) 
shows that it is possible for <j> to increase at the cen- 
ter of the star, and then to turn around and decrease at 
larger r. This turnaround can occur because the pressure 
decreases with r, causing the force to increase. 

Once the star is large enough that (f> c ~ (f>t , chameleon 
"screening" makes the field far less responsive to further 
increases in GM+/r+. Screening is the stellar analog of 
the chameleon thin shell effect, in which the deviation of 
<f> from its background value is sourced only by a thin 
shell of matter near the surface of an object. The thin 
shell effect becomes important when the field at the cen- 
ter of an object approaches the minimum of its effective 
potential inside that object, in a precise analogy to stellar 
screening. For a sufficiently large star, the effective stel- 
lar mass that acts as a source to (f> will be much smaller 
than the actual stellar mass. 

Specifically, we define the relativistic analogues of the 
bare (linear) and effective (screened) masses: 



-^screen 



47rr 2 (p- 3P)dr, 



dr 



3r 2 dV cS 

o 2G dej) 

2/ , 3r 2 dU 

W(,-3P) + -- 



(16) 



dr. (17) 



In the linear regime, when the slope of the potential is 
small, M screen ?s Mu near . Since these are cases where 
P/p <C 1, Miinear ~ A/* and the field feels the total mass, 



6 




FIG. 4: Chameleon screening for two values of x\, with n = 1 and p c = lOOpA- Left: The chameleon field's deviation <j> c — (f>as 
from its background value rises with gravitational potential until GM*/r* w 1 — 0as, after which the chameleon becomes 
significantly less responsive to further increases in potential. The thin curve denotes negative values of (f> c — 4>dS, reached at 
very high pressures, P c > p c . Right: As the chameleon enters the non-linear regime, its source M screen decreases with respect 
to its linear regime source Mii near . This screening becomes important around GM*/r* w 1 — 0ds- 



which is also the source of the gravitational potential 
GM+fr*. Therefore the change in the field is proportional 
to the gravitational potential (see Sec. IIV Al for a more 
extended treatment). 

As the potential and the change in the field become 
large, the field rolls to a steeper part of the potential. 
Thereafter the field source is screened by the potential, 
reducing the source from Mi; near to M screen . As C — ^ 0t, 
this screening becomes complete at the center of the star, 
and the chameleon is only sensitive to stellar matter at 
larger r. Fig.|3]illustrates screening in two different ways. 
As the potential GM+/r+ is increased in Fig. [4] (left), 
4>c ~ 0dS increases steadily until GM+/r± »1- 4>dS- Be- 
yond that point, <fi is insensitive to further increases in the 
potential. Fig. 0] (right) shows that the onset of this in- 
sensitivity coincides with the decrease of M screen /Mii near . 
As cj) c — > 4>t, the chameleon "sees" a smaller and smaller 
portion of the star, so further increases in the potential 
are unable to push the field all the way to the curvature 
singularity. 

There has been much confusion in the literature be- 
tween non-linear chameleon effects and relativistic ef- 
fects. As Fig. @] makes clear, screening is a chameleon 
effect that is totally unrelated to strong gravity. The 
chameleon enters the non-linear regime when GM+/r* m 
4>t ~~ 0ds ~ 1 — 4>dS and screening becomes important. 
This is true even when 1 — 0ds «C 1, for which a star with 
GM^/r^ w I — 0ds is non-rclativistic. fn the Starobin- 
sky f(R) model, I — 0ds is determined by the choice of 
model parameters n and x\, for n = 1, 1 — 0ds = X\ 2 ■ 
Chameleon effects and relativistic effects will coincide 
when x\ is of order unity, as in much of the literature. 
For x\ ^S> 1, chameleon effects will appear in objects 
with potentials much smaller than unity, fn fact, the 
appearance of chameleon effects for the galactic poten- 
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FIG. 5: Potential and chameleon field in high-pressure stars. 
As P c is increased beyond ^p c , the potential reaches a maxi- 
mum and (f> c — 4>dS becomes negative. 



tial is required for solar system tests of gravity [9(. fn 
other words only Starobinsky models with ij > 1 and 
1 _ 0ds ^ 10~ 6 — 1CT 5 are in fact viable. 

Finally, we note that Fig. Q] (left) hints at a truly rel- 
ativistic effect. For x\ = 3.6, the plot of C — 0ds begins 
to decrease as GAh/r* is increased beyond about 0.3, 
corresponding to P c > p c /3. Note that general relativity 
would require negative scalar curvature R at the center 
of such a star. For stars of this size, Mi; near begins to 
decrease with respect to M+. As we continue to increase 
Pc/ Pc beyond 1/3, we find that Mii near and M scrccn can 
become negative. This means that <f> c — </>ds will be neg- 
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ative in a sufficiently high-pressure star, as in the thin 
curve in Fig. 0] (left). Furthermore, we find that GM+/r+ 
reaches a peak in the high P c regime, and then begins to 
decrease with P c , as shown in Fig. [5] 

The high central pressures P c /p c 3> 1 do not actu- 
ally reflect the average P/p c through the star. In fact, 
for these stars, the stellar radius, as well as the pressure 
and field profile near the surface, become independent of 
P c /p c , keeping GM+/r* constant. The rapid rise in the 
central pressure in a small central core is accompanied 
by a suppression of the central field <fi c — ► which al- 
lows potential and pressure gradients to be balanced in 
hydrostatic equilibrium. 

For n = 1 and x\ = 3.6, we find that the gravita- 
tional potential reaches a maximum of GM*/r* = 0.345 
at P c /p c = 0.78. This is in contrast to general relativity, 
in which GM^/r^ increases monotonically toward 4/9 as 
Pel pc — ► oo. (For the stars considered here, corrections 
to this general relativistic upper limit due to the presence 
of a de Sitter horizon are negligible [13] ■) However, such 
high pressures are not thought to be realized in any star 
composed of ordinary matter, so this difference between 
general relativity and f(R) gravity is unlikely to be useful 
for observationally distinguishing between these theories. 



IV. ANALYTIC ARGUMENTS 

In the previous section, we used numerical examples 
for several choices of f(R) parameter values in order to 
show that: 

1. the chameleon field (f> can be in one of two regimes, 
linear and non-linear; 

2. the linear regime is characterized by a linear scaling 
between the field and the gravitational potential, 
0c - 0ds = GM*/r*; 

3. the non-linear regime is characterized by a field pro- 
file that increases near the center, but turns around 
before hitting the singularity; 

4. the transition between these regimes occurs when 
the gravitational potential approaches 1 — 0dS i the 
distance in field space between the de Sitter value 
and the singularity, and is unrelated to relativistic 
effects. 

Here, we use analytic arguments to generalize these 
four results to a broader range of parameter values and 
stellar densities. In Sec. IIV A[ the linear field equation 
is solved exactly for a non-relativisitic, constant density 
star in order to show that <p c — 0ds — GM+fr*. The 
linear approximation breaks down rapidly as GM^/r^ 
approaches 1 — 0dSj even in a non-relativisitic star. In 
Sec. IIVB1 we study the field near the stellar center in the 
non-linear regime. 4>{r) increases slowly in response to 
the pressure, even as —GM(r)/r decreases. The field ap- 
proaches the singularity, but turns around before reach- 
ing it. We also use our analytic solution to show that 




0.077 

1e-06 1e-05 0.0001 0.001 0.01 
r [Gpc] 

FIG. 6: Numerical (</> n um) and approximate (<^ un scr) field pro- 
files for a star in the unscreened limit. We assume n = 1, 
x\ = 3.6, p c = 100pA, and P c /p c = 10 -4 . The dotted and 
solid vertical lines correspond to the stellar surface and the 
horizon, respectively. 



matching the exterior boundary condition requires expo- 
nential precision in setting the central field value. The 
implications of this exponential tuning are studied for 
different densities and parameter values, such as those 
used in llfjl. 



A. Unscreened stars 

Throughout this section, we work in the non- 
relativistic limit, P c /p c , GM+fr* <C 1, in which the equa- 
tion of motion (|5|) reduces to 



,., 2 , dU 8ttG 
r dip 6 

Approximating the potential by 

dU 2 t± ± \ 
— w -m dS (<p- 0ds), 



valid for <fi near 4>&S, wc find 



[r) - (pds = 



3m dS 

/~t e 



- <->int " 



ih(m d sr) 



r < rv, 
r > r*. 



(18) 



(19) 



(20) 



The constants Cj nt and C ext arc found by matching 
4>(r) and 4>'(r) at the stellar surface, r = r*, and the 
resulting field profile is shown in Fig. [SI At the stellar 
center, 



4>c - 0dS 



GM* 



2(1 - (1 +?7i ds r,)e- md ^.) 



("ids?"*) 2 



(21) 



The quantity in brackets approaches 1 in the limit of 
small mdsr*, which is an excellent approximation because 
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"Ms ~ Ha- Thus we have 4> c — 0ds = GM+fr* for an 
unscreened non-relativistic star. Since GM+/r± <C 1, this 
also validates the assumption in Eq. (TiU]) . 



B. Screened stars 

As GM+/r* approaches 1 — fa§, 4>c —> F so (TT51) implies 
that becomes large. Thus, our approximation that 
ui(p ~ mas is violated; the unscreened scaling breaks, 
and the field no longer responds linearly with the po- 
tential. In this screening limit, the full field profile can 
only be calculated numerically. On the other hand, in 
the stellar interior, where screening is nearly complete, 
analytic solutions are available that provide insight into 
the numerics. 

To zeroth order, the field near the stellar center sits at 
the minimum of its effective potential, 



d 



V eS (r, 4> min ) 



dU 



8ttG 



(p c -3P) = 0. (22) 



This is because the field at r w is insensitive to the field 
outside the star, which is many Compton wavelengths 
away, and can adjust itself to minimize the local effective 
potential. Note that </> m i n (0) = fa; fa n - m is a generaliza- 
tion of fa to r > 0. In particular, we see that 4> m in{r) 
increases with r at the center toward the singularity, be- 
cause P(r) decreases. Since the field must eventually 
match onto the exterior solution, with cf> (r*) < 0, the 
field must turn around. This applies to non-relativistic 
as well as relativistic stars. Given a fixed pressure profile 
it is straightforward to solve for m i n - In the high cur- 
vature, R 3> Rq limit, the minimum corresponds to the 
general relativistic expectation that R = 8irG(p — 3P), 
so analytic expressions for P(r) are available. 

One can solve the field equation © iteratively to ob- 
tain successively better approximations to the screened 
solution. To first order the screened solution fa cr be- 
comes 



dU 



isr_ 

2N 



-) 

2B J 



D, 

(23) 

where N and B are given by the general relativistic so- 
lution. 

We stop at first order, since the true solution must 
depart substantially from the screened solution in the 
outer regions of the star in order to match the exterior 
solution smoothly. Fig. [7] shows fa cr compared with the 
numerical solution of the equations of motion, for a star 
with P c /pc = 0.1 in an f(R) model with n = 1 and 
x\ = 3.6. Notice that the two solutions only deviate in a 
shell of mass near the stellar radius, corresponding to a 
region where the integrand of M screen becomes nonzero. 
Thus the outer regions of the star source the deviation of 
from fa cr that allows it to roll continuously to fas ■ This 
numerical solution simply reflects a smooth interpolation 




0.0001 



FIG. 7: 1 — <j> for the numerical computation, the screening 
limit c^scr, and perturbations around it 4> acl + 8(j> in the con- 
stant mass limit of Eq. (|28|l . We assume n = 1, xi = 3.6, 
p c = 100pA, and P c /p c = 0.1. The dotted and solid verti- 
cal lines correspond to the stellar surface and the horizon, 
respectively. 



between the interior screened solution and the exterior 
boundary condition. 

Deviations from the screened solution at the outskirts 
imply that even deep in the interior there will be small 
deviations. These deviations 5<j) = cf> — fa cl - are governed 
by 



N' 






4> SCI +8<t> 










2N ' 




d(j) 





(24) 



where the N and B solutions can be iterated to the ap- 
propriate order. 

Note that in the approximation that the field has rolled 
only a small distance from its central value, \5<f>\ <C 1— fa, 

dU 4>scr+> 



where 



d<j)/dR 



m t 5<fi 



R 



(25) 



(26) 



Equation (f2"4")l therefore becomes source free and has 
Yukawa-like solutions which exponentially grow with 
nitr. Since mtr* 3> 1, this implies that <p c must be ex- 
ponentially close to, but not exactly equal, fa cr at r = 0, 
in order for the deviations from the screening solution 
to become significant only near the stellar radius. Thus 
numerical solutions for stars that are screened for a sub- 
stantial part of their interior are difficult to find numeri- 
cally by shooting from a central value fa. 

These considerations can be made more concrete for 
non-relativisitic stars. In this case, the pressure pro- 
file P(r) w P c (l — r 2 /r+ 2 ), and the corresponding ze- 
roth order solution m i n contains a quadratic piece. The 
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screened solution deviates from d m in as 



(0scr 0min) — 0min "i" 0min ~ 30 m i n 

r 



2m t 2 



(27) 



Furthermore, deviations away from the screened solution 
grow as 



S(f>(r) = 8cj){Q) 



sinh(m t r) 



m t r 



(28) 



There are a number of interesting properties of the nu- 
merical solution that can be gleaned from this analytic 
treatment. We have already shown that (f> t = </>min(0) 
is the threshold between field solutions that decrease 
monotonically (4> c < (f>t) and those that increase at 
low r (</> c > <f>t). Now we see from Eq. (|28| that, if 
4> c > </> S cr(0) = 0t + k, then the field will monotonically 
increase. Thus there is an interval of width k in field 
space for which the field will increase at the center and 
then turn around. For a star with n = 1, x\ = 3.6, 
p c = 100/9A, and P c = 0.1p C7 the width of this interval is 
of order 10 -9 . Fig. [7] also shows </> scr + 5<fi for this case 
under the approximation of Eq. (|28p . with 6(f)(0) chosen 
so as to match the turnaround point from the numerical 
solution. Note that deviations from this approximation 
are expected after the field turns around to match the 
exterior boundary since the constant mass approxima- 
tion breaks down. We also drop the small relativistic 
correction here for simplicity. 

Furthermore, most of this interval 4>t < 4>c < <^t + t 
corresponds to fields that turn around at r ~ m%~ , much 
earlier than in a typical screened star. If we want this 
turnaround position rt ur n, at which (rturn) = 0, to be 
much larger, then much more tuning in the field value 
is necessary. For r turn > ro, the lower bound on 5(f)(0) 
becomes 



50(0) > ~-K(m t r ) 2 e- m < r <>. 



(29) 



We see from our numerical solution Fig. [3] (right) that 
rturn ~ r*/2 for P c /pc = 0.1. From (|29[) . the interval in 
field space for which r tU m > r*/2 is of order 10 -39 , ap- 
proximately the amount of tuning that we needed earlier 
in order to find a numerical solution. We can estimate 
the numerical difficulty of finding (f> c though the shooting 
method, for a general star, by expressing k in terms of 
the f(R) model parameters n and x\, 



9-2 8,l+3 (2n+l) 2 

-(4n+2) 



1 + (n + l)x\ 



(1 + x\) n+1 



PA 



(30) 



At p c = 2 x 10 6 pa, a star with the same m t r* as the one 
discussed above, and rturn > r*/2, will require that (f> c be 



tuned to a precision 10 65 , and this tuning only becomes 
worse at higher n. 

V. CONCLUSIONS 



We have studied static, spherically symmetric stars in 
f(R) theories of modified gravity. Such theories behave 
like ordinary gravity with a chameleon field, a matter- 
coupled scalar with non-linear self interactions. We have 
found the chameleon field profile inside a star by numer- 
ically integrating the modified Einstein equation. The 
critical features of these numerical solutions are exposed 
by analytic arguments. 

We find solutions to the equations of motion over a 
wide range of central pressures, 10 -5 < P c /p c < 1000, 
and gravitational potentials, 1.4 x 10~ 5 < GM+/r+ < 
0.345 in Starobinsky's model of f(R) gravity. This range 
of potentials extends from small, sun-like stars all the way 
to highly relativistic neutron stars. Thus the existence of 
relativistic stars cannot be used to rule out f(R) theories 
of gravity. 

We have shown that, in the non- linear regime of the 
chameleon field, stars are screened. That is, (f> is sourced 
only by the outer portion of a star's mass, analogous to 
the chameleon thin shell effect. This keeps (f> from reach- 
ing the singularity s ing = 1 as GM^/r^ is increased. As 
the stellar density increases at fixed GM*/r+, the dispar- 
ity between the Compton wavelength of the field and the 
stellar radius makes numerical solutions unfeasible, but 
does not imply that solutions do not exist. 

Screening sets in when GM+/r± ps 1 — ^jSj that is, 
when the gravitational potential approaches the distance 
in field space between the singularity and the de Sitter 
background. This marks the transition between the lin- 
ear and non-linear regimes of the chameleon. This tran- 
sition is unrelated to strong gravity. Its appearance at 
GM 4 /r t f=a 10 _1 in KM [l(| is an artifact of their choice 
of an f(R) model with 1 — 0ds ~ 10 _1 . We have explored 
models where 1 — (f>ds ~ 10 -3 and confirmed that screen- 
ing sets in at GM t /r t r* 10~ 3 , that is, in non-relativistic 
stars. 

The fact that chameleon effects can appear at much 
lower potentials than strong gravity is essential to the 
construction of viable f(R) theories @. Viability re- 
quires that the Galaxy, with potential $ ~ 10~ 6 — 10~ 5 , 
be in the non-linear chameleon regime. It is only in this 
regime that deviations from general relativity in the solar 
system are sufficiently suppressed to satisfy local tests. 
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